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We report anomalous quantum transport features in bilayer graphene in presence of a random distribution of 
structural vacancies. By using an efficient real-space Kubo-Greenwood transport methodology, the impact of a 
varying density of dimer versus non-dimer vacancies is investigated in very large scale disordered models. While 
non-dimer vacancies are shown to induce localization regimes, dimer vacancies result in an unexpected ballistic 
regime whose energy window surprisingly enlarges with increasing impurity density. Such counterintuitive 
phenomenon is explained by the formation of an effective linear dispersion in the bilayer bandstructure, which 
roots in the symmetry breaking effects driven by dimer vacancies, and provides a novel realization of Dirac 
semimetals in high dimension. 
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Introduction.- Single layer graphene (SLG) has attracted 
a great attention owing to its remarkable electrical, chem¬ 
ical and mechanical properties providing an endless list of 
novel opportunities for practical applications 111]. SLG pos¬ 
sesses a linear electronic spectrum with chiral (A-B sublat¬ 
tice) symmetry which lead to exotic low-energy trai^ort 


such as Klein tunneling QS), weak antilocalization |U,l5t], 
half-integer quantum Hall effect Hi. 

Bilayer graphene (BLG) differs from SLG by the 
parabolic band dispersion which however retains the chi¬ 
ral nature of low-energy electronic excitations. One of the 
salient and unique property of BLG is the possibility of cre¬ 
ating an electronic bandgap by applying external gate volt¬ 
age lili . However surprisingly, the understanding of quan¬ 
tum transport in disordered BLG remains far less understood 
than the SLG case, because of the enhanced structural com¬ 
plexity. Experimental studies evidence critical differences 
in transport behaviors of SLG and BLG US- Some gener¬ 
alization of the localization theory in BLG has been derived 
d , while a minimum conductivity a min — jh is pre¬ 
dicted at the charge neutrality point (CNP) Il2l - [fl . A recent 
scanning tunneling microscopy (STM) study shows that va¬ 
cancies in graphite induce peculiar impurity states known as 
zero-energy modes (ZEMs) which are maximally localized 
at the defect position and then decay as the inverse of the 
distance from the vacancy d. The impact of ZEMs on 
BLG is so far poorly understood, especially the role played 
by dimer and non-dimer vacancies which strongly differ in 
terms of symmetry breaking characteristics. 

In this Letter, we start by analysing the localization fea¬ 
tures of ZEMs in BLG for all types of vacancies, and found 
a highly inhomogenous sublattice state population (pseu¬ 
dospin polarization) as reported in STM experiments 115 1. 
The depletion of low energy states in one sublattice is addi¬ 
tionally further classified into two different classes depend¬ 
ing on the vacancy position. Then by using efficient com¬ 
putational methods, we explore quantum transport in disor¬ 


dered BLG and analyze how the local nature of the vacancy 
(dimer versus non-dimer) impact on scattering and localiza¬ 
tion phenomena. The type and concentration of vacancies 
are found to dictate the nature of the transport regime which 
ranges from weak localization to anomalous ballistic con¬ 
duction. 

Electronic features induced by vacancies in graphene.- 
BLG (EiglHa) and (c)) can be considered as two coupled 
SLGs with the top layer (in red) shifted a carbon bond from 
the bottom layer (in black). Consequently, BLG consists of 
four carbon atoms in its unit cell, two carbons Ai, Bi in 
the bottom SLG unit cell and A 2 , B 2 in the top layer where 
B 2 lies on the top of Ai, namely dimer or a sites whereas 
Bi,A 2 are called non-dimer or /3 sites. The tight-binding 
Hamiltonian model for BLG reads ll^E3] 


H =-7o + 71 X + 73 X (1) 

(ij),I—1,2 i i,i' 

- 74 X{“m“ 2 ./ + h.c. + {a+.ai^i + bl^b 2 ,^} 

i,j' i 


where ( = 1,2 labels the bottom and top layer respectively. 
The annihilation (creation) operators acting on Ai, Bi, A 2 , 
B 2 are denoted ai, bi, 02 , 62 (oi^, a^, b^). The first 

term in TL describes the intralayer hopping between nearest- 
neighbor TT-orbitals. The second term denotes the interlayer 
hopping (71 = 340 meV) between dimer sites {Ai,B 2 }, 
while the third term gives the interlayer coupling between 
Bi and its closet A 2 with 73 = 280 meV. The fourth term 
corresponds to the hopping from Ai to its nearest A 2 site 
and from Bi to its nearest B 2 neighbors (74 = 145 meV). 
The energy asymmetry between dimer and non-dimer sites 
is taken into account by introducing in the final term of TL 
an on-site energy difference A = 9.6 meV between dimer 
sites. All these parameters are derived from the ab-initio 
calculations |3 1 ^ . 


The Local Density of States (LDOS) of pristine BLG are 
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first scrutinized (dash curves in the inset of Figlljb) or (d)). 
Since the two layers are identical one can restrict the discus¬ 
sion to a (dimer) and /3 (non-dimer) sites. The LDOSs on 
both sites show a sudden change in the slope si E = ± 71 , 
especially for a sites (black dashed line) owing to the con¬ 
tribution of higher energy bands ll^ . The LDOSs also 
clearly exhibit fingerprints of pseudospin polarization on 
each layer, in the sense that the state mainly populate /3 sites 
(LDOS^g ^ 0 and LDOS„ = 0). 




FIG. 1: (color online): Sketch of a dimer (a) and non-dimer (c) 
vacancy (open dashed circle) in a BLG composed of a top (in red) 
and a bottom (in black) layers. State distribution around a dimer (b) 
and a non-dimer (d) vacancy on its layer in the opposite sublattice. 
Power law fit (red solid curve). Inset: LDOSs shown in (a) and (c), 
closest to dimer (b) and non-dimer vacancies (d), together with the 
pristine BLG (dashed curves). 

A monovacancy on BLG is simulated by removing a car¬ 
bon atom from the pristine BLG which can be either chosen 
on the dimer or the non-dimer site. FigUla) and (c) show the 
sketch of BLG with dimer and non-dimer vacancies (dashed 
circle), while the decay of the impurity states are given in 
Fig- Eb) and (d), respectively. We first consider the ef¬ 
fect of a single vacancy on its layer (here the bottom layer). 
Our results (FigEb) and (d) inset) show that both dimer and 
non-dimer vacancies have a strong impact on the layer where 
the vacancy resides with a strong depletion of ZEM on the 
vacancy sublattice and a state abundance in the other sub¬ 
lattice (similarly to SLG, see Supp. Mat. lEoll ). This is seen 
in FiglT](b) and (d) (insets) by comparing the LDOSs of the 
two nearest sites in the bottom layer (black and red solid 
lines) with the pristine BLG case (dashed lines). Dimer va¬ 
cancies also more strongly impact at high energy as seen in 
the decay of the LDOS at E = ±71. The spatial distribution 
of the sublattice abundant state induced by the vacancy in 
its layer is further investigated in FigEb) and (d) where the 
data (closed circle) for both dimer (FiglT](b)) and non-dimer 


vacancy (FiglT](d)) are fitted by power law (red lines). 
The depletion of the charge density in the vacancy sublattice 
can be actually classified into two different classes. The first 
class involves six second-nearest neighbors of the vacancy 
and forms a hexagonal lattice (green dashed line in FigEc)) 
with an enlarged graphene lattice constant by a factor of v/3, 
whereas the second class together with the vacancy populate 
the centers of the previous hexagons. Such long range nature 
of the impurity state distribution will be key in understand¬ 
ing anomalous ballistic transport. We next investigate the 
effect of a single BLG vacancy on its adjacent layer. While 
the dimer vacancy introduces ZEMs on both sublattices of 
the adjacent layer (the peaks at i? = 0 in green and blue 
curves in EigEb) inset), the non-dimer vacancy leaves the 
second layer almost unaffected. Indeed, the LDOSs for a 2 
(blue solid line) and j32 (green solid line) site in EigEd) in¬ 
set are almost identical to the ones for pristine BLG (dashed 
lines). 

Transport properties in the dilute vacancy limit.-Change 
transport properties of BLG are investigated for a finite den¬ 
sity of vacancies, differentiating the cases with only dimer or 
non-dimer vacancies (uncompensated cases) from the case 
with equally distributed mixture of both types (compensated 
case). Here we also assume an equal distribution among top 
and bottom layers. We use a real-space order-N wave packet 
evolution approach ll^ . The Kubo-Greenwood conduc¬ 
tivity is written as a{E,t) = e^p{E)AX^{E,t)/t, where 
p{E) is the DOS and AX^{E, t) is the mean quadratic dis¬ 
placement of the wave packet and gives the diffusion co¬ 
efficient D{E,t) = AX'^{E,t)/t. With disorder, D(t) 
changes from a ballistic motion to a saturation regime, from 
which the mean free path £e is deduced through te{E) = 
D'^™[E)/2v{E) (v{E) is the velocity and the maxi¬ 
mum value). At long times, D{E, t) eventually decay owing 
to quantum interferences which drive the system either to 
weak or strong (Anderson) localization regime m. 

EigEa) shows the total DOS for BLG with 0.05% of 
dimer vacancies a (red solid curve), non-dimer (3 (blue solid 
curve) and mixed a — jd case (black solid curve), together 
with the total DOS for pristine case (dashed curve). As al¬ 
ready suggested in Eig[T] (b) and (d) (insets), especially at 
the CNF, dimer vacancies more strongly affect the DOS than 
non-dimer ones, whereas the depletion of state at E = ±71 
is unobservable. However, a marked difference is observed 
on the mean free path 4 (FigEc)) and semiclassical conduc¬ 
tivities Use = e^p{E)D'^^ (EigEd)) with two significant 
peaks at E = ±71. The values of 4 (resp. a sc) at such en¬ 
ergy are higher for the dimer than for the non-dimer vacan¬ 
cies by a factor of 3 (resp. a factor of 2) indicating a stronger 
scattering efficiency for non-dimer vacancies. This likely 
originates from the imbalance of state at E = ±71 on the 
two sublattices around the non-dimer vacancies (EigEb), (d) 
insets). The stronger impurity scattering strength at high en¬ 
ergy for non-dimer vacancies is also confirmed by the fact 
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FIG. 2: (color online)(a) Total DOS for BLG samples with 
0.05% vacancies together with the pristine DOS (dashed curve), 
with only dimer vacancies (red curve), non-dimer vacancies 
(blue curve) or a mixed case (black curve), (b), (c) and (d) 
D{E = 0,f), Ir{E) and asc{E) corresponding to the DOSs 
in (a), (d, inset) crsc(0.02%)/(Tsc(0.05%) (black curve) and 

o'sc(0.05%)/crsc(0.5%) (red curve). 


that le and age are almost the same for the compensated 
case (black solid curves in Figl2c),(d)) and the non-dimer 
vacancy case (blue solid curves). For three cases are al¬ 
most the same but age for the dimer vacancy is larger due to 
enhanced DOS. D{E^ t) further evidence much stronger lo¬ 
calization effects for the dimer vacancies (FiglJJb)), consis¬ 
tently with the DOS results (FiglT](b) and (d)). Importantly 
in the low impurity density limit (< 0.05%), regardless the 
nature of vacancies, a localization regime is obtained for the 
whole energy spectrum. 

However, surprising features are obtained for high enough 
vacancy density when the affected spatial areas around va¬ 
cancies start to overlap. One first scrutinizes the ratios of 
semiclassical conductivity age for the compensated case in 
the dilute limit with 0.02%, 0.05% of vacancy, and dense 
limit with 0.5% of vacancies (FiglJJd) (inset)). In the dilute 
limit age perfectly obeys the Fermi’s golden rule ~ l/rii 
in the whole energy band (black solid curve in FiglJJd) 
(inset). The Fermi’s golden rule, however, underestimates 
(overestimates) the values of age at CNF {E = ± 71 ) in the 
dense impurity limit (red solid curve in FiglJJd) (inset), i.e. 
when increasing the vacancy coverage age decreases faster 
at i? = ± 7 i and slower at the CNR 

Transport properties in dense vacancy limit. -FigI3 and 
FiglD show results for 0.5% of non-dimer vacancies and 
compensated case (Fig[3 ) and dimer vacancies (Fig®. 
Charge transport in presence of non-dimer vacancies and for 
the compensated case similarly exhibit a localization behav¬ 
ior (as seen in Fig[3a)). The low-energy semiclassical con- 


FIG. 3: (color online) D{t) (a) and asc (b) for 0.5% vacancy in 
the non-dimer (red curves) or compensated case (black curves), 
with the minimum conductivity Ae^/h for BLG (dashed curve), 
(c) Act(L) = a{L) — CTsc for 0.5% of non-dimer vacancy (open 
circle) and the fit to weak localization theory (red solid curve). 


ductivities also clearly saturates at ~ 4e^//i (dashed 

curve) il4 l. 

The calculation of the quantum correction of the semiclas¬ 
sical conductivity Aa{L) = a(L) — age at CNP (Figl^c)) 
further confirms that the non-dimer vacancies induce local¬ 
ization, since a very good fit Aa{L) ~ a^ln(L/(*) is 
obtained with a = —0.295 ^ — I/tt in agreement with 2D 


weak localization theory (A(t(L) ~ —^ln(L/(*)). 

The case of high density of dimer vacancy is remarkably 
different from the non-dimer one. Fig®a) shows D{E, t) 
for 0.5% dimer vacancies at CNP (red solid line) and at high 
energy E = O.O 570 (red dashed line). While high energy 
charge transport quickly enters the localization regime (see 
D{E = 0 . 0570 , t)), in sharp contrast D{E = 0, t) increases 
linearly without any sign of saturation, as expected in a truly 
ballistic motion. This remarkable behavior is captured in the 
quantum conductivity scaling a{t) in the vicinity of CNP. 
Fig®b) shows that a{t) at f = 5 ps (dashed lines) is smaller 
than its value at f = 40 ps (solid lines) around CNP. A 
crossover from the ballistic regime to the localization regime 
actually occurs at a critical energy Ee = O.OO 870 (red dotted 
lines). This anomalous behavior is also observed for a higher 
dimer vacancy density of 5% (Fig|4ja, b) (blue lines))., for 
which a significant enlargement of the ballistic-transport en¬ 
ergy range is obtained (Ee = O.OI 770 (blue dotted lines)). 
This is a remarkably counterintuitive phenomenon in which 
the quantum conductivity increases with disorder density. 
As explained below, the interlayer hopping between non¬ 
dimer sites 73 governs such anomalous behavior. First, by 
setting 73 = 0 for the case of 5% dimer vacancies, D{t) is 
seen to exhibit a strong localization with a very fast decay at 
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FIG. 4: (color online) (a) D{E = 0, t) (solid lines) and D{E — 
0 . 0570 , t) (dashed lines) for 0.5% (red lines) and 5% (blue lines) 
of dimer vacancies. The green line gives the case of 5% of dimer 
vacancy with 73 = 0. (b) Conductivities at f = 5ps (dashed lines) 
and at t = 40ps (solid lines) for 0.5% (red lines) and 5% (blue 
lines) of dimer vacancies, (c) Charge distribution on the vacancy 
sublattice around a dimer vacancy, which separate into two classes 
(see text), together with their fittings to the power law (solid lines). 


long times and a value reduced by two orders of magnitude 
when compared to the case with 73 0 (Figlla) (green 

solid line)). 

To rationalize such effect, let us then consider a simple 
tight-binding model by suppressing the terms 74 and A in 
Eq.(l 2 ]i, which is equivalent to the continuum model 22 1 


H = -- 


2 m 


0 

0 


+ 


0 TT 
TT^ 0 


( 2 ) 


It describes two possible interlayer hopping events between 
non-dimer sites, the ones with dominant charge density 
(FigHJb) and (d) inset (dashed lines)). The first term acti¬ 
vates hopping via dimer sites which generate a mass term 
TO = 7 i/2u^ and a parabolic band i^. The second term de¬ 
scribes a direct hopping 73 between non-dimer sites in the 
bottom and top layers with velocity ?;3 = {'\/Z/ 2 )ay^/Ti. 
This term is dominant at low energy and generates a linear 
energy dispersion Sin 


. = ^(rapP-^cos(3« + (|i) (3) 

where ^ = 4-1 (— 1 ) corresponds to the K {K') valley and 
(j) labels the momentum direction. This linear dispersion oc¬ 
curs for e < 571(^3/?;)^. Hereafter we will show that this 
value is enhanced in presence of dimer vacancies and yield 
anomalous ballistic motion. 


The impurity state created by a dimer vacancy is mainly 
localized on non-dimer sites, and has almost no weight on 
dimer sites around the defect. Such sublattice polarization 
induced around the vacancy is shown in FiglHb). Fig| 4 ](c) 
further evidences the depletion of FDOS on the dimer sites 
of the layer where the vacancy is lying. Similarly to the SFG 
case (see 1 ^ ). the FDOSs on dimer sites can also be sepa¬ 
rated into two classes (circle and square symbols) which are 
both well fitted with the scaling law p = pQ — br~^-^ (black 
solid lines) where po is the FDOS far away from the vacancy. 
From Fq.(l 2 ]i, one observes that the abundance of charge den¬ 
sity on the non-dimer sites leads to the enhancement of the 
second term in Fq.© which corresponds to the direct hop¬ 
ping between non-dimer sites. The reduction of electronic 
state on the dimer sites, on the other hand, leads to a reduc¬ 
tion of the first term in Fq.(| 2 ]i. Indeed, this term involves a 
hopping process in which electron has to hop into a dimer 
site. This process is limited by dimer vacancies because it 
creates around it, in dimer sites, a wide area of depletion 
of electron density, i.e. prohibiting electron residence on 
the dimer sites in the vacancy vicinity. We can thus con¬ 
sider BFG with a high enough density of dimer vacancies 
as an effective pristine BFG with renormalized parameters 
V* < V and > V3. The more dimer vacancies in the BFG, 
the smaller the renormalized value of v* and the larger the 
renormalized value of ^3. This thus leads to the expansion 
of the energy range (571 (^3/^*)^ ) where the linear disper¬ 
sion dominates. This also explains the increase of critical 
energy Ec with the dimer vacancy density in FigHfb). Such 
chiral electrons in the renormalized linear BFG energy dis¬ 
persion share many similarities with low-energy excitations 
propagating in SFG (the sublattice index is replaced by the 
layer index) yielding backscattering suppression. This situ¬ 
ation thus provides a novel realization of a Dirac semimetal 
in high dimensionality, which is an unconventional transport 
regime provoked by symmetry effects of impurities. This 
could motivate a systematic exploration of quantum trans¬ 
port in irradiated BFG, following pioneering experimental 
studies il 5 l. 
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